Lattice Boltzmann method at finite-Knudsen numbers 
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A modified lattice Boltzmann model with a stochastic relaxation mechanism mimicking "virtual" 
collisions between free-streaming particles and solid walls is introduced. This modified scheme 
permits to compute plane channel fiows in satisfactory agreement with analytical results over a 
broad spectrum of Knudsen numbers, ranging from the hydrodynamic regime, all the way to quasi- 
free flow regimes up to Kn ~ 30. 

PACS numbers: 



The dynamic behaviour of flows far from hydro- 
dynamic equilibrium is an important subject of non- 
equilibrium thermodynamics, with many applications in 
science and engineering. The non-hydrodynamic regime 
is characterized by strong departures from local equi- 
librium which are hardly handled on analytical means. 
Consequently, much work is being devoted to the devel- 
opment of computational techniques capable of dealing 
with the aforementioned non-perturbative and non-local 
effects. Recently, the lattice Boltzmann (LB) method 
has attracted considerable interest as an alternative to 
the discretization of the Navier-Stokes equations for the 
numerical simulation of a variety of complex flows 0. 
Extending LB methods to non-hydrodynamic regimes is 
a conceptual challenge on its own, with a variety of mi- 
crofluidic applications, such as flows in micro and nano 
electro-mechanical devices (NEMS, MEMS) 0. De- 
partures from local equilibrium are measured by the 
Knudsen number, namely the ratio of molecular mean 
free path, to the shortest hydrodynamic scale, Ih'- 
Kn = Im/lh- Ordinary fluids feature Kn < 0.01, while 
high Knudsen numbers are typically associated to rar- 
efied gas dynamics, where Im is large because density is 
small, typical case being aero-astronautics applications. 
More recently, however, the finite-Kn regime is becoming 
more and more relevant for a variety of microfluidics ap- 
plications in which the Knudsen number is large because 
of the increasingly smaller size of the devices. It is also 
worth emphasizing that the finite-Knudsen regime may 
also bear relevance to the problem of modeling fluid tur- 
bulence [2b i^J . Recent work indicates that LB may offer 
quantitatively correct information also in the finite-Kn 
regime Isf. This hints at the possibility that LB may 
complement or even replace expensive microscopic sim- 
ulation techniques such as kinetic Monte Carlo and/or 
molecular dynamics. On the other hand, this is also 
puzzling because finite-Knudsen fiows are in principle 
expected to receive substantial contributions from high- 
order kinetic moments, whose dynamics is arguably not 
quantitatively correct because of lack of symmetry of the 
discrete lattice 0, 0- In this work we point out that, 
irrespectively of symmetry considerations, any LB ex- 
tension aiming at describing non-hydrodynamic flows in 



FIG. 1: A simplified sketch showing the velocity structure 
for the D2Q9 model in a plane channel flow. In LBE mod- 
els an hypothetic particle forming an angle 6 with respect to 
the X axis (population "A") is counted amongst the particles 
in population "1". This particle would then, in the absence 
of collisions (i.e. for high- Knudsen numbers), be forced to 
stream parallel to the x axis without ever impacting on the 
walls. In real flows instead, its natural fate would be to move 
along a straight line until a collision with the wall would oc- 
curr at position "W". 



the finite-Knudsen regime must first provide an adequate 
description of the dynamics of free-streamers, namely 
ballistic particles whose motion, in the limit of infinite 
Knudsen number, escapes any thermalization by either 
bulk or wall collisions. 

This problem is connected to a pathology of the LB 
model, and in fact, of any discrete velocity model rep- 
resenting a whole set of molecular speeds within a finite 
solid angle, by a single discrete speed. The practical con- 
sequence of this drastic reduction of degrees of freedom 
in momentum space is a very efficient numerical scheme, 
but also the pathology connected to the fact that parti- 
cles moving parallel to the wall (for example in a channel 
flow) never impact on the boundaries, so that under the 
effect of pressure gradients or body forces within the fluid 
flow, they may enter an unrealistic free-acceleration (run- 
away) regime. More generally, runaway must be expected 
whenever the particle mean free path exceeds the longest 
free-flight distance allowed by the geometrical set up. 

In this work we propose a mechanism to obviate 
this runaway pathology. Our method is demonstrated 
through the numerical simulations of a force-driven plane 
Poiseuille flow, as compared with existing analytical the- 
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ories 

The lattice Boltzmann method has been described in 
detail in many publications, and here we shall only re- 
mind a few basic ideas behind this technique. The sim- 
plest lattice Boltzmann equation looks as follows 0: 

f,{x + Atc,,t + At)-f,{x,t)= (1) 

where fi{x, t) = f{x, v = Ci,t), i = 1, ... ,n, is the prob- 
ability of finding a particle at lattice site x at time t, 
moving along the lattice direction defined by the dis- 
crete speed Ci, and At is the time step. The left-hand 
side of this equation represents free-streaming, while the 
right-hand side describes collisions via a simple relax- 
ation towards local equilibrium /f (a local Maxwellian 
expanded to second order in the fluid speed) in a time 
lapse T — Lu~^. This relaxation time fixes the fluid kine- 
matic viscosity as = c^(t — 1/2), where Cg is the sound- 
speed of the lattice fluid, {cg = l/-\/3 in the present 
work). Finally, Fi represents the effects of an external 
force. The set of discrete speeds must be chosen such 
that mass, momentum and energy conservation are ful- 
filled, together with rotational symmetry. Once these 
symmetries are secured, the fluid density p — fiy ^^'^ 
speed u = ^^fiCi/p can be shown to evolve according 
to the (quasi-incompressible) Navier-Stokes equations of 
fluid-dynamics. In this paper, we shall refer to the nine- 
speed D2Q9 model shown in Figure 

While at low Knudsen local equilibrium is established 
via bulk collisions, in the regime Kn > 1, the leading 
equilibration mechanism is provided by the interaction 
with the boundaries. Various types of boundary condi- 
tions have been used for LB simulations of finite-Knudsen 
flows 0, 0, 0, 01 . In this work, we confine our atten- 
tion to bounce back (BB) Q and kinetic, as recently 
introduced by Ansumali and Karlin (AK) [l3|. The for- 
mer implement no-slip flow speed at the wall via particle 
reflections, while the latter are based on the idea of rein- 
jecting particles from the wall according to a local equi- 
librium distribution with wall speed and temperature. 
While AK were explicitly designed with finite-Knudsen 
effects in mind, this is clearly not the case for BB. At 
high-Knudsen, neither BB nor AK are expected to pro- 
vide realistic results, since neither of them has been de- 
signed to do so. 

It is therefore of interest to compare these two quite 
distinct types of boundary conditions throughout the full 
range of Knudsen numbers. 

In a intermediate/high-Knudsen flow the main mech- 
anism for thermalisation is provided by collisions with 
solid walls. This mechanism is non-local, meaning by this 
that, in order to be thermalized, a particle residing at a 
given lattice site (see population "A" in Fig. needs to 
reach the wall (position "W" in Fig. where it is scat- 
tered along some different direction, thereby transferring 



some momentum to the wall. From Figure ^ it is clear 
that in the LB representation, all molecules with momen- 
tum in the solid angle [—tt/& : 7r/8] around population 
"1", are collapsed into population "1". In the limit of 
infinite Knudsen numbers, population "1", being parallel 
to the walls, behaves like a ballistic beam, whose motion 
escapes both bulk collisions and wall scattering. Due to 
the presence of an external drive, these beam particles 
enter a free-acceleration regime, which ultimately leads 
to runaway behaviour in time. Therefore, in the way to 
a finite-Knudsen LB scheme, no matter how accurate in 
terms of tensorial symmetry, a mechanism to thermalize 
these ballistic streamers needs to be introduced. 

Here we propose a mechanism to mimick non-local wall 
collisions through virtual wall collisions (VWC for short). 
At every time step each free-streamer is assigned a virtual 
speed (called "A" in Fig. {vx,Vy) = c-(cos 9, sin 6*), with 
uniformly distributed in [— 7r/8, 7r/8]. Assuming colli- 
sions are distributed according to a Poissonian, we define 
p as the probability to "hit" (i.e. to collide at least once) 
the wall in a single time step, multiplied by the proba- 
bility that no molecular collisions occurred during that 
time step. In practice, this means that a wall-collision 
event is performed at each time step and each lattice site 
with probability: 

p(.,,;0-e-V^"-(l^e- ''"°-T-» ) (2) 

where the first factor is the probability of undergoing no 
bulk collision during the flight of (average) length H to 
the wall, when the average mean-free-path for bulk col- 
lisions is A (i.e. —1/Kn ~ —X/H). The second term 
describes the fact that after a flight of length L the par- 
ticle hits the wall with probability one. This is mod- 
eled at each time step in the following way. In a time 
step dt the particle moves a distance dl — cdt; since 
on average there is one collision every L lattice spac- 
ings, the probability of a wall collision in a time step dt 
is given by exp (—dl/L) — exp(— cdi sm(6)/H) (in our 
units cdt = 1). 

Note that p goes to zero in the limit of Kn 0, so 
that virtual wall collisions fade away in the hydrody- 
namic regime, as they should. Due to the non-analytic 
dependence in Kn the hydrodynamic limit is recovered 
faster than any power of Kn. 

The probability p goes to zero also in the limit 9^0, 
corresponding to the case of the virtual speed falling back 
into a free-streamer (see Figure 1). However, by design, 
the occurrence of such a limit has now a vanishingly small 
probability, much like in a continuum flow. 

Upon hitting the wall, the free-streamer feeds the span- 
wise moving populations according to the rules adopted 
for true wall-collisions (see Figure 1). For instance, for a 
virtual collision of a right-moving particle with the top 
wall, we have: 

/( = /l(l-p), /7.8 = /7,8+p/l/6, /l = /4 + 2p/i/3 
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where the coefficients {1/6,4/6,1/6} correspond to the 
weights of the nine-speed local equilibrium. 

Thanks to this mechanism, momentum is removed 
from populations aligned with the wall boundaries and 
re-distributed along the direction orthogonal to the 
boundaries. This allows the flow to relax towards a (non- 
local) equilibrium also in the limit of high Kn numbers. 

We simulate a two-dimensional flow in a rectangular 
duct of size L along the streamwise direction (x) and 
width H < L along y. The flow is driven by a volumet- 
ric force (force per unit volume) F^. = pc^jj^ along the 
streamwise direction. In the hydrodynamic limit, such 
force leads to a parabolic velocity profile of amplitude 
Uq. Both BB and AK boundary conditions are applied 
on-site. The Knudsen number is tuned by changing the 
value of the viscosity, according to the expression: 

Kn = l^/H = v/{c,H) 

Knudsen numbers in the range 10~'^ < Kn < 30 have 
been simulated, at a fixed Mach number Ma — u/cs — 

0. 03. The grid resolution used was 101 x 21, although a 
few simulations at higher resolution of 101 x 41 have also 
been performed to check the grid-independence of our re- 
sults. After discarding the initial transient time, velocity 
profiles were averaged in time in order to minimize statis- 
tical fluctuations associated with the stochastic nature of 
the VWC mechanism.. The simulation was run until the 
relative change of the profiles fell below a given threshold 
(comparable with machine precision). 

One of the major successes in the early days of kinetic 
theory was the prediction of a minimum of the mass flow 
rate as a function of the Knudsen number around Kn ^ 

1. Interestingly, both BB and AK boundary conditions 
exhibit a minimum mass flow rate around Kn ~ 1 (see 
Figures [3 and . This indicates that, contrary to the 
common tenet, LB can capture the so-called "Knudsen 
paradox" . However, they both severely overestimate the 
mass flow in the high-Kn region, Kn > 1. This is due to 
the presence of a finite fraction of ballistic streamers, as 
discussed earlier on. Figures [21 and 13 show the behaviour 
of the normalized mass flow rate Q = Ux{y) / {'ov), for 
the case with and without virtual wall collisions. 

The numerical results are compared with the analytical 
expression obtained by Cercignani in the zero and infinite 
Knudsen limits respectively (both fiuxes are normalized 
in units of Qpv): 

Qo = l/{6Kn) + s + {2s^ ~1) Kn (3) 
^ 7r-i/2 log (Kn) (4) 

where s = 1.015 (see |1]). 

By increasing the Kn number the flow slips more and 
more, until the centerline speed becomes comparable 
with the sound speed, so that the LBE simulation breaks 
down. 
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FIG. 2: AK boundary contions: normalized mass flux as a 
function of Kn number, with (o) and without (t) virtual 
wall collisions. The dashed line is Cercignani's prediction ^ 
for low-Kn. Dotted lines are prediction Q for high-Kn mul- 
tiplied by factors 1,2 and 4, from bottom to top. Following 
Cercignani, the mass flow is normalized with 6 pi'. 
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FIG. 3: BB boundary contions; normalized mass flux as a 
function of Kn number, with (o) and without (t) virtual 
wall collisions. Notations are the same as in the Figure |5| 



As shown in Figure |21 by introducing virtual wall col- 
lisions, a nearly quantitative agreement with analytical 
results is observed. It is worth stressing that no free pa- 
rameter has been used in the simulations. The highest 
achievable Kn number may not be high enough to test 
the asymptotic prediction. Therefore, in order to appre- 
ciate the quality of our results, we have plotted Cercig- 
nani's prediction, Qcxj, in Figures |2] and [SJ together with 
its magnifications by factors 2 and 4. 

Back in 1867, Maxwell predicted that the slip length 
of a rarefied gas flowing on a semi-infinite slab should 
be proportional to the mean free path. Is — l.lAQKn, 
[l7| |. Nearly four decades ago, Cercignani computed the 
next order corrections in Kn, leading to the following 
quadratic expression for the slip speed: 

Vs = aKn + bKn^ (5) 
where a = 1.146 ■ (¥] and 6 = -0.9075 • 

The dependence of the slip speed on the Knudsen num- 
ber is presented in Figure^ The slip velocity is defined 
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FIG. 4: Wall slip velocity as a function of the Knudsen num- 
ber for AK (plusses) and BB (triangles) boundary conditions. 
The solid line is the wall slip velocity extrapolated from the 
values of the flux computed assuming a flat velocity distribu- 
tion across the channel. The dashed line is Cercignani's low- 
Knudsen analytical prediction lO . Error bars as described in 
the text. 



situations, such as non ideal geometries, high shear rates 
[l^ and thermal effects jljf . Another interesting point to 
be explored for the future is the potential benefit of using 
multi- relaxation [2(1 El| and entropic LB schemes for 
a better description of fluid- wall interactions. 

The authors thank S. Ansumali and I. Karlin for dis- 
cussions and for making available their implementation of 
AK boundary conditions in a early stage of this project. 
SS acknowledge discussions with A. Garcia. Finan- 
cial support from the NATO Collaborative Link Grant 
PST.CLG.976357 and CNR Grant (CNRCOOBCBF-001) 
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as the value (at the wall position) of the parabolic fit of 
the averaged velocity profiles. Boundary data were ex- 
cluded from the fit. The error bars on the slip velocity 
were estimated as the discrepancy between the results of 
the 21 X 101 and 41 x 101 simulations. 

From Figurc2]it is clearly seen that BB boundary con- 
ditions fail to reproduce the slip velocity. Thus, our re- 
sults support previous criticism on the use of BB bound- 
ary conditions at finite-Knudsen 0,^^. However, ki- 
netic boundary conditions provide quantitative agree- 
ment with analytical results up to Kn 1, even beyond 
the theoretical expectations. At intermediate and high 
Knudsen both BB and AK provide similarly inaccurate 
results. This is because the physics becomes more and 
more insensitive to the BB or AK boundary conditions 
and dependent on the VWC mechanism instead. Indeed, 
as virtual collisions are introduced, both BB and AK pro- 
vide satisfactory agremeent with analytical data. The 
solid line in Figure 01 corresponds to the ratio Q/{pH) 
and the fact that this quantity matches almost exactly 
the measured slip speed indicates that the flow profile is 
basically flat, as predicted by the asymptotic theory. Fi- 
nally, our results do not extend below Kn = 0.01 simply 
because an accurate evaluation of the very small values 
of Vs in this regime requires a much higher transverse 
resolution than adopted here. 

Our results indicate that a standard nine-speed LB 
scheme equipped with Ansumali-Karlin boundary condi- 
tions and a virtual wall collision mechanism, can capture 
salient features of channel flow in both hydrodynamic 
and strongly non-hydrodynamic regimes. Of course, this 
is only the first step towards a systematic inclusion of 
high-Knudsen effects in the lattice kinetic framework, 
and much further work is needed to address more general 
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